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Chapter  1 

INTRODUCTION 


This  report  describes  the  details  of  an  algorithm  which  allows  to  find  an  iso-surface 
for  a  given  three-dimensional  property  distribution.  Input  to  the  algorithm  may  be  a 
set  of  scalar  values  given  at  regularly  spaced  grid  points.  Output  will  be  triangulated 
polygonal  approximation  of  the  iso-surface  for  a  specified  property  value. 
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Chapter  2 

THE  CONCEPT  OF 
ISO-SURFACES 


The  concept  of  iso-surfaces  assumes  that  a  scalar  property  p  (for  example:  pres¬ 
sure,  temperature,  or  magnitude  of  velocity)  is  distributed  inside  a  three-dimensional 
volume  where  it  is  continuously  defined  as  p  =  p(x,  y,  z).  For  our  considerations,  we 
assume  that  the  volume  of  interest  is  a  cuboid  (a  right  parallelepiped)  and  bounded 
by  rectangular  faces.  Therefore,  the  property  distribution 

p  =  p(*>  y,*) 

is  defined  for  the  given  ranges  of  x,  y,  and  z  : 


®mtn  —  %  5; 

Vmin  —  V  <  Jfmui 
*mi»  ^  *  $  ^mas  • 

Inside  the  volume,  the  scalar  property  p  changes  continuously  from  point  to  point 
and  exhibits  values  within  some  problem  dependent  range 
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Pmin  —  P  —  Pmax • 


To  study  the  character  of  the  three-dimensional  property  distribution,  a  constant 
property  value  pc  can  be  chosen  and  all  points  exhibiting  this  property  value  pc  can 
be  found  inside  the  given  volume.  The  points  of  constant  property  value  pc  "an  be 
determined  from 

P(z,y,z)  ~Pc  =  0, 

which  is  an  implicitly  defined  three-dimensional  surface.  Such  a  surface  is  called  an 
iso-surface  and  represents  the  locus  of  all  points  where  the  property  p  assumes  the 
given  constant  property  value  pe.  The  geometry  of  an  iso-surface  is  problem  dependent 
and  varies  from  very  simple  shapes  (for  example:  a  plane,  a  sphere,  or  a  cylinder)  to 
extremely  complex  forms.  Sample  iso-surfaces  are  illustrated  in  Figure  2.1.  As  can 
be  seen  from  Figure  2.1,  a  single  iso-surface  may  consists  of  two  or  more  unconnected 
surface  parts. 

Note  that  within  a  given  property  distribution,  an  infinite  number  of  iso-surfaces 
can  be  found,  each  of  them  characterized  by  a  different  constant  property  value  pc. 
Figures  2.1a  and  2.1b  illustrate  two  iso-surfaces  for  two  different  pe  values  within  the 
same  volume. 

For  a  given  pe,  an  iso-surface  exists  inside  a  volume  as  long  as 

Pmin  —  Pc  ^  Pmax- 

If  pe  is  outside  the  range  of  Pmin  -  Pmax,  then  an  iso-surface  for  the  given  pic  cannot 
be  found  within  the  given  volume. 
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Chapter  3 

PRINCIPLES  OF  ALGORITHM 


3.1  Numerical  Definition  of  Property  Distribu¬ 
tion 

The  algorithm  which  is  described  in  the  following  requires  a  numerical  definition  of 
the  property  distribution  p  =  p(x,  y,  z)  inside  a  given  volume.  Discrete  property  values 
Pij,k  must  be  provided  at  the  intersection  points  of  a  rectangular,  three-dimensional 
mesh  of  straight  lines.  Figure  3.1  illustrates  a  volume  with  a  sample  regular  grid. 
The  grid  is  defined  by  the  arrays 

Xi ,  i  from  1  to  NX, 

Yj ,  j  from  1  to  NY, 

Zky  k  from  1  to  NZ. 

Note  that  this  definition  allows  for  variable  spacing  between  the  grid  lines  in  all  three 
directions.  The  intersections  of  the  grid  lines  are  the  grid  points  defined  as  the  vectors 

GP  ij*  =  (XitYjtZk). 

The  total  number  of  grid  points  NG  for  a  given  volume  is 
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xl  x2  x3  x4 


Figure  3.1:  Computational  volume  with  Grid  lines 


NG  =  NX  x  NY  x  NZ. 

At  each  of  the  grid  points  GPjj  *,  the  scalar  property  value  P%jjt  is  given.  The  range 
of  the  property  values  for  a  given  property  distribution  can  be  found  from 

PMIN  =  MIN  (PiJtkt  i  =  l,NXJ  =  l,NY-,k=l,NZ) 

Pm  ax  =  MAX  {Pijjkt  i  =  l,NX,j  =  ltNY-k  =  1,  NZ) 

3.2  Subdividing  the  Volume  into  Boxes 

The  regular  grid  can  be  used  to  divide  the  given  volume  into  sub- volumes.  Each 
sub-volume  is  again  a  cuboid  or  right  parallelepiped  and,  for  simplicity,  will  be  called 
box  or  grid  box  in  the  following.  A  box  extends  between  two  consecutive  grid  lines  in 
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Figure  3.2:  Local  Identification  System  of  a  Box 

all  three  directions,  i.e.,  a  box  is  bounded  by  the  six  planes  X  =  X{  and  X  =  Xl+j, 
Y  =  Yj  and  y  =  ,  and  Z  —  Zk  and  Z  =  Zi+i  .  The  total  number  of  boxes  NB 

inside  the  given  volume  is 

NB  =  (NX  -  1)  x  (NY  - 1)  x  (NZ  -  1). 

Figure  3.2  illustrates  single  box  and  introduces  a  local  identification  system  for 
the  geometric  elements  of  a  box: 

•  A  box  is  defined  by  the  8  vertices  V*,  k  =  1,8.  The  vertices  coincide  with  the 
given  grid  points  GP,jt*.  At  each  vertex  V*,  the  property  value  Pu  is  known. 

•  The  8  vertices  form  the  12  edges  Ei ,  l  =  1,12. 

•  The  12  edges  define  6  rectangular  faces  Fm,  m  =  1,6. 
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Figure  3.3:  Iso-surface  patch  within  a  Box 

3.3  Box-by-Box  Surface  Construction 

To  find  the  iso^surface  for  a  given  constant  property  value  Pc,  a  box- by-box  surface 
construction  method  is  used.  Each  of  the  NB  boxes  is  investigated  individually.  If 
the  iso-surface  intersects  with  single  box,  a  patch  of  the  iso-surface  which  lies  inside 
the  box  can  be  found  as  illustrated  in  Figure  3.3.  First,  the  intersection  points  S, 
(i  =  l,  NS)  of  the  iso-surface  with  the  edges  of  the  box  are  found  assuming  a  linear 
distribution  of  the  property  value  along  each  edge.  Next,  the  intersection  points  S< 
are  connected  by  line  segments.  Each  line  segment  lies  in  one  of  the  faces  of  the  box. 
The  line  segments  build  the  boundaries  of  the  iso-surface  patch  inside  the  box  and 
define  a  closed  polygon  which,  in  general,  is  non-planar.  As  shown  later,  the  number 
of  intersection  points  NS  may  vary  between  3  and  12.  Therefore,  the  shape  of  the 
resulting  polygons  may  range  from  a  simple  triangle  to  a  12-sided  polygon.  In  some 
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Figure  3.4:  Two  boxes  sharing  &  face 


cases,  more  than  one  polygon  may  result  from  the  process.  All  resulting  polygons 
can  be  decomposed  into  triangles.  These  triangles  are  used  as  the  basic  geometric 
elements  for  the  iso-surface  representation. 

The  entire  iso-surface  is  constructed  by  combining  the  surface  patches  found  within 
each  of  the  NB  boxes  of  the  volume.  If  the  iso-surface  extends  from  one  box  into 
a  neighboring  box,  the  transition  is  always  continuous:  Both  boxes  share  a  common 
face  and  the  polygons  from  both  boxes  share  the  same  edge  lying  in  that  face.  This 
is  illustrated  in  Figure  3.4.  Therefore,  the  sum  of  all  polygons  (or  the  sum  of  all 
triangles  decomposed  from  the  polygons)  represent  the  entire  iso-surface  inside  the 
given  volume. 

Not  all  of  the  NB  boxes  of  the  volume  contribute  to  the  construction  of  the 


surface.  An  individual  box  is  not  intersected  by  the  iso-surface  if 


Pc  >  ft.  *  =  1,8 

or  if 

*  =  1,8. 

Such  a  box  is  an  empty  box.  Moreover,  the  entire  given  volume  may  be  empty  if 


Pc  >  Pm  ax 


or 

Pc  <  Pmin- 

In  this  case,  no  iso-surface  exists  inside  the  entire  volume  for  the  given  value  Pc. 


3.4  Polygonal  Approximation  of  Iso-surface  Patch 
within  a  Box 

Three  levels  of  approximations  are  applied  to  the  construction  of  iso-surface  patch 
inside  each  grid  box.  First,  the  intersection  points  S,  are  found  using  linear  inter¬ 
polation  along  each  edge.  Second,  boundary  curves  of  the  iso-surface  patch  are  ap¬ 
proximated  by  straight  lines.  Finally,  the  iso-surface  patch  is  approximated  by  sum 
of  triangles. 

The  bilinear  interpolation  can  be  applied  on  a  face  of  a  box  to  show  a  part  of 
above  approximation  process.  This  interpolation  can  be  used  to  find  out  traces  of  the 
iso-surface  on  the  face. 
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Bilinear  interpolation  is  defined  by  four  corner  points[3] 


P(u,t>)  :  P(0,0),  P(0,1),  P(1,0),  and  P(l,l) 

The  corner  points  define  linear  boundary  curves  (straight  lines)  for  the  patch  which 
is  formed  by  the  interpolation.  Any  points  on  the  patch  is  linearly  interpolated  by 

Q(u,  v)  =  P(0,0)(1  -  u)(l  -  v)  +  P(0,1)(1  -  u)v  +  P(1,0)«(1  -  v)  +  P(l,l)uv 


For  one  face  of  a  box  we  can  use  the  more  familiar  x,y  and  z  coordinates  where  z 
coordinate  represents  the  property  value  P  at  a  point  defined  by  x  and  y.  Also,  the 
bilinear  interpolation  can  be  rewritten  in  matrix  form 


Q(*.»)  = !  (l  -  x)z  1 


P(0, 0)  P(0,1) 
P(1.0)  P(l,l) 


(1-v) 

s 


or 


Q(^.y)  =  [  (1  -x)(l  -y)  (l-x)y  *(l-y)  xy  ] 


P{0,0) 

P(0,1) 

P(1.0) 

P(U) 


Figure  3.5  shows  one  bilinear  patch  defined  by  four  corner  points 


P(0,0)  =  (0,0,5), 

P(0, 1)  =  (0,1,3), 

P(1,0)  =  (1,0,2)  and 
P(M)  =  (1.1.7) 

We  can  observe  that  all  the  iso-parametric  lines  (  x  =  const,  or  y  =  const)  are 
straight. 
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Figiire  3.5:  A  bilinear  surface  patch 

The  contour  line  can  be  defined  by  Qt(xty)  s  ft  where  P  -  tk  • 

»  made  the  range  of  P  th„  •  '  “  U  Pc 

5  oi  s,  there  is  a  contour  i;n.  a. 


line  is 


contour  line,  and  the  equation  for  the 


contour 


c  =  *(0,0)  -  l(*(0,0)-*(I,0)).y(/,(00)  /,(oi)) 

+  *^,(0,0) 

We  rewrite  above  equation  with  substitutions  of  foUowings 

P’M  =  ft* 

^O.l)  =  ft,, 

^•c.0)  =  *n  and 

^0.1)  =  *j 


-  ft.,)  -  r(PM  _  ft,  _  +i>  ), 

,jW  ~  F«»  ~  *(flw  -  *,)  -  ft 
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We  rearrange  it 


_ (Pojo  ~  Pc)  +  x(-Pq,q  4-  Pi, o) _ 

^  (Po.o  —  Po.i)  +  x(~Pq,o  +  Po,i  +  P 1,0  —  P 1,1) 


In  short  form 


where 


and 


If 


a  +  bx 
c  +  dx 


a  =  P o.o  —  Pc 

b  =  —Po,o  +  Pifl 

c  =  Po.O  —  Po,l 

d  =  —  Pofl  +  Pp,l  +  P\fl  —  P\,X 


c  +  dx  ^  0 


c  +  dx  -  0, 


then 


a 

y  =  oo  at  x  =  -- 


As  an  example,  let’s  investigate  a  bilinear  patch  shown  on  Figure  3.5  where  Poja  — 
5,  Po, i  =  3,  Pifi  =  2,  and  P14  =  7.  We  want  to  find  out  contour  line(s)  of  property 
value  of  4.  Then, 

a  =  1 

b  =  -3 
c  =  2 

d  =  -7 
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Figure  3.6:  Contour  lines  of  a  Bilinear  surface  patch 


The  contour  line  is  defined  by 

1  -3x 
V  ~  2  -  7x 

The  contour  line  is  shown  as  Figure  3.6.  We  can  observe  that  the  contour  line  is  not 
straight  whereas  the  iso-parametric  lines  are  all  straight. 

Another  example  case  which  shows  the  process  of  polygonal  approximation  is 
generated  which  is  defined  by  a  function 


P  =  0.04xy  +  z 


Figure  3.7:  Real  iso-surface  patch 
within  a  box  defined  by  following  eight  vertices 

(-50,-50,-50)  where  P  =  +50 
(+50,  —50,  —50)  where  P  =  —150 
(—50,  +50,— 50)  where  P  =  —150 
(+50,  +50,  —50)  where  P  =  +50 
(—50,  —50,  +50)  where  P  =  +150 
(+50,  —50,  +50)  where  P  =  —50 
(—50,  +50,  +50)  where  P  =  —50 
(+50,  +50,  +50)  where  P  =  +150 

Figure  3.7  shows  the  box  and  an  iso-surface  of  value  Pc  =  0  within  the  box.  The 
iso-surface  will  be  approximated  to  a  sum  of  triangles  by  the  algorithm  developed 
in  this  report.  First,  the  intersection  points  of  the  iso-surface  with  edges  of  the  box 
will  be  approximated  by  linear  interpolation  along  each  edge.  On  the  edge  Ex  (see 
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Figure  3.8:  Approximated  iso- surface  patch 


Figure  3.2),  y  =  -50  and  z  =  —50.  The  exact  *  value  of  the  intersection  point  is 
calculated  by 


0  =  0.04x(— 50)  +  (-50) 


and 


x  =  —25 

The  linear  interpolation  gives  x  value  of  the  intersection  point  as 

x  =  -25 


The  top  face  and  bottom  face  have  four  intersection  points  each.  By  connecting  these 
intersecting  points  with  straight  lines,  we  approximate  the  iso-surface  patch  within 
the  box  as  figure  3.8.  This  non-planar  8-sided  polygon  will  be  divided  into  triangles 


Chapter  4 

CLASSIFICATION  OF 
BOX- SURFACE 
INTERSECTIONS 


As  can  be  seen  from  Figures  3.3  and  3.8,  an  iso-surface  may  intersect  a  selected 
box  in  various  ways.  Before  an  algorithm  can  be  devised  to  find  the  boundaries  of 
the  surface  patches  (in  the  form  of  non-planar  polygons),  the  possible  topologies  of 
the  polygons  must  be  investigated  individually. 

4.1  The  Total  Number  of  Intersection  Cases 

The  points  S  of  the  polygons  are  found  by  intersecting  the  iso-surface  with  the 
edges  of  the  box.  Assuming  a  linear  change  of  the  property  value  along  each  edge, 
the  intersection  point  is  found  from 

S  =  V,  +  (VJ-V1)x^^  (4.1) 

where,  Vi  and  V2  denote  the  vertices  at  either  end  of  a  selected  edge,  P\  and  P7  are 
the  corresponding  properties  given  at  Vj  and  V2,  respectively.  Figure  4.1  illustrates 
the  finding  of  S  along  the  edge. 
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Figure  4.1:  Linear  interpolation  along  an  edge 

The  different  patterns  of  box-surface  intersection  are  determined  by 
of  edge  intersections  found  and  by  the  topological  configuration  of  the 
points.  An  edge  may  have  no  intersection  point  at  all  if 

Pc  >  Pi  and  Pc  >  Pj 

or 

Pc  <  Pi  and  Pc  <  Pj 

An  edge  may  have  one  intersection  point  somewhere  in  the  middle,  may  have  a  sin¬ 
gle  intersection  point  at  either  of  both  limiting  vertices,  or  the  entire  edge  may  lie 
completely  on  the  iso-surface  if  Pc  =  Pi  =  Pi- 

The  various  edge  intersections  are  created  by  different  combinations  of  property 
values  at  the  eight  vertices.  As  obvious  from  the  above,  only  three  conditions  at  each 
vertex  are  significant:  the  property  at  a  vertex  is  either  grater  than  Pc,  is  less  than 


the  number 
intersection 
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Pc,  or  is  identical  with  Pc-  Therefore,  a  total  of 

3*  =  6,561 

cases  are  possible  for  a  single  box.  Because  of  symmetry,  many  of  these  6,561  cases 
can  be  handled  as  identical  patterns.  Nevertheless,  the  number  of  remaining  cases  is 
too  large  for  the  development  of  an  effective  algorithm. 

4.2  Reduced  Number  of  Intersection  Cases 

The  possible  number  of  box-surface  intersection  can  be  signi.icantly  reduced  by 
eliminating  the  condition  that  a  property  at  a  vertex  is  identical  with  Pc .  In  practical 
applications,  the  probability  for  an  equality  condition  is  extremely  small.  In  the  rare 
case  that  such  an  equality  occurs,  a  vertex  property  which  is  identical  with  Pc  can 
be  modified  (for  example,  lowered)  by  a  very  small  amount.  The  algorithm  uses  the 
following  method  for  a  consistent  property  modification: 

PTOL  ■  (Pm ax  -  Pmin)  *  0.00001 

c  10  *  min.  positive  floating  point  value  with  4  bytes 

rlmin  ■  l.E-37 

IF  (PTOL  .LT.  rlmin)  PTOL-  rlmin 
c 

IF  (  ABS(Pl-PC)  .LT.  PTOL)  PI  -  PC  -  0.5*PT0L 

Figure  4.2  shows  a  result  of  such  modification  for  a  vertex  with  four  neighboring  faces. 
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Figure  4.2:  Slight  modification  of  property  value 

After  this  modification,  only  two  significant  conditions  remain  for  a  vertex:  the 
property  value  is  either  greater  than  Pc  or  is  less  than  Pc-  Therefore,  a  total  of 

2*  =  256 

intersection  cases  remain  and  need  to  be  investigated.  It  also  can  be  considered  as 
bit  patterns  of  8  bits.  Each  bit  represents  each  vertex  of  a  box  and  it  can  have  two 
values,  0  and  1,  depending  on  the  property  values  there. 

4.3  Classification  in  Groups 

A  study  of  all  the  256  cases  of  box-surface  intersections  shows  that  the  number 
of  edge  intersections  can  be  used  to  classify  the  cases  into  nine  different  groups. 
Table  4.1  shows  those  9  groups  with  their  frequencies  and  unique  patterns.  The 
minimum  number  of  edge  intersections  is  zero,  and  the  maximum  is  twelve.  Note 
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Figure  4.3:  Sample  cases  of  Nine  groups 

that  cases  of  1,  2,  10,  and  11  edge  intersections  do  not  exist.  Figure  4.3  shows  sample 
cases  of  nine  groups.  All  the  256  cases  grouped  by  the  number  of  edge  intersections 
are  enumerated  in  Appendix  A. 

In  the  following,  possible  connecting  patterns  for  each  group  are  identified.  One 
pattern  is  identified  as  unique  by  its  topology,  i.e.  number  of  its  sides  (edges)  and 
number  of  separate  parts  within  a  box.  Figure  4.4  shows  those  unique  patterns  of 
each  group. 

•  Group  1 :  No.  of  edge  intersections  =  0 ;  2  cases  of  the  256  cases  have  no  edge  in¬ 
tersections.  All  values  of  8  vertices  are  higher  than  the  iso-surface  value  or  lower 
than  that.  One  case  is  the  dual  case  of  the  other.  The  bit  pattern  of  00000000 
produces  the  same  topological  pattern  with  the  bit  pattern  of  11111111.  In  these 
cases,  no  surface  patches  or  polygons  exists,  but  it  is  counted  as  one  pattern. 
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(see  the  following  page  28.1) 


Figure  4.4:  Unique  patterns  of  iso-surface  patches  of  9  groups 

•  Group  2  :  No.  of  edge  intersections  =  3  ;  16  cases  of  the  256  cases  have  3  edge 
intersections.  In  these  cases,  there  is  only  one  unique  pattern.  The  pattern  is 
formed  by  connecting  3  edge  intersections  in  any  order.  They  form  a  triangle. 

•  Group  3  :  No.  of  edge  intersections  =  4  ;  30  cases  have  4  edge  intersections. 
There  is  one  unique  pattern.  Four  edge  intersections  form  a  4-sided  polygon. 

•  Group  4  :  No.  of  edge  intersections  =  5  ;  48  cases  have  5  edge  intersections. 
There  is  only  one  pattern.  Five  edge  intersections  form  one  5-edged  polygon. 

•  Group  5  :  No.  of  edge  intersections  =  6  ;  64  cases  have  6  edge  intersections. 
There  are  two  different  patterns.  One  of  them  is  6-edged  single  polygon.  The 
other  pattern  is  made  of  two  triangles. 

•  Group  6  :  No.  of  edge  intersections  =  7  ;  48  cases  have  7  edge  intersections. 
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Figure  4.4:  Unique  patters  of  iso-surface  patches  of  9  groups 


28.1 


There  are  two  patterns.  One  of  them  is  7-edged  single  polygon.  The  other  has 
one  triangle  and  one  quadrilateral. 

•  Group  7  :  No.  of  edge  intersections  =  8  ;  30  cases  have  8  edge  intersections. 
There  are  three  different  patterns.  One  of  them  are  8-edged  single  polygon. 
Another  pattern  has  one  triangle  and  one  5-edged  polygon.  Still  another  has 
two  quadrilaterals. 

•  Group  8  :  No.  of  edge  intersections  =  9  ;  16  cases  has  9  edge  intersections. 
There  are  three  patterns.  One  of  them  is  composed  of  9-edged  single  polygon. 
Another  has  one  triangle  and  a  6-edged  polygon.  The  last  pattern  is  composed 
of  three  triangles. 

•  Group  9  :  No.  of  edge  intersections  =  12  ;  2  cases  have  12  edge  intersections. 
There  are  four  patterns.  1)  Four  triangles.  2)  Two  triangles  and  one  6-edged 
polygon.  3)  Two  6-edged  polygons.  4)  One  12-edged  polygon. 
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Chapter  5 


HANDLING  OF  AMBIGUOUS 
CASES 


If  any  of  6  faces  of  a  grid  box  has  4  edge  intersections,  there  is  no  unique  way  of 
connecting  edge  intersections  on  the  face.  This  produces  ambiguities  on  the  algorithm. 
One  must  find  the  correct  way  of  connecting  intersection  points.  For  a  face,  there  are 
3  possible  ways  of  connection.  Figure  5.1  shows  a  face  with  ambiguity  and  possible 
connections. 

One  simple  solution  is  to  subdivide  the  domain  (i.e.  &  box)  into  smaller  cells  until 
the  ambiguity  is  dissolved.  This  is  called  as  adaptive  subdivision ,  or  tesselation.  But 
this  method  requires  functional  values  at  internal  points  of  one  grid  box,  which  are  not 
directly  availabe  in  a  given  discrete  data  set.  That  requires  interpolation  of  proper 
form  (local  or  global).  Also,  it  creates  cracks  on  an  iso-surface  between  different 
size  cells,  i.e.,  between  subdivided  box  and  undivided  box.  Because  straight  line  is 
used  to  connect  the  intersection  points,  the  side  of  the  iso-surface  patch  found  on  the 
smaller  size  box  does  not  match  with  that  on  the  bigger  size  box.  The  subdivided 
box  has  more  intersection  points  and  they  approximate  the  iso-surface  more  closely. 
Figure  5.2  shows  one  face  shared  by  two  boxes,  one  of  them  is  subdivided  into  eight 
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Figure  5.1:  A  face  with  Ambiguity 

smaller  boxes. 

Triangulation  of  the  grid  box  is  another  solution.  A  triangular  face  can  have  ex¬ 
actly  two  intersections,  if  any.  There  is  only  one  way  of  connecting  them.  Figure  5.3 
shows  a  triangular  face.  In  general,  no  ambiguities  exists  if  a  simplex  is  the  partition¬ 
ing  cell.  A  simplex  is  the  simplest  linear  decomposition  of  n-space.  For  3-D  space,  it 
is  tetrahedron.  Figure  5.4  shows  possible  patterns  within  a  tetrahedron.  This  method 
requires  triangulation  of  domain  (into  tetrahedra),  and  also  requires  function  values 
at  number  of  internal  points.  Finding  intersection  points  along  edges  of  triangulated 
cell  is  compuationally  more  expensive  than  along  edges  of  box  shaped  cells,  where 
edges  can  be  aligned  with  coordinate  axes. 

Still  another  solution  is  utilizing  function  values  at  the  center  of  ambiguous  faces. 
From  the  fact  that  the  iso-surface  divides  the  higher  value  region  from  the  lower  value 
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Figure  5.4:  Iso-surface  patterns  within  Tetrahedra 

region,  the  correct  way  of  connecting  points  can  be  determined.  The  function  value 
at  the  center  of  ambiguous  face  is  tested  against  the  iso-value  Pc .  If  the  centroidal 
value  is  higher  than  Pc,  the  two  vertices  with  lower  values  will  be  isolated,  vice 
versa.  Figure  5.5  shows  a  face  with  four  edge  intersections  and  the  centroidal  value. 
This  method  gives  the  same  result  as  when  the  box  is  subdivided  into  24  tetrahedra, 
if  additional  information  is  only  used  to  get  the  correct  connecting  information,  i.e., 
additional  intersections  created  inside  and  on  the  box  are  neglected.  This  method  still 
needs  function  value  at  the  centroid  of  the  ambiguous  face.  Some  form  of  interpolation 
should  be  applied. 

Linear  interpolation  (specifically,  mean  of  values  of  4  vertices)  is  used  in  the  im¬ 
plementation  to  get  the  center  value.  It  is  a  simplified  form  of  the  distance  weighted 
interpolation.  The  example  case  shown  in  Figure  3.8  is  constructed  by  using  the 
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Figure  5.5:  Four  edge  intersections  and  the  Cetroid  value 

mean  value  to  get  the  connecting  information  on  top  and  bottom  faces.  If  the  mean 
of  property  values  of  four  vertices  is  used,  the  centroid  value  of  the  top  face  is  higher 
than  the  iso-value  in  Figure  3.8.  That  of  the  bottom  face  is  lower  than  the  iso-value. 

As  linear  interpolation  is  also  used  for  finding  edge  intersections,  it  is  consistent 
throughout  the  domain.  It  will  produce  consitent  results  at  neighboring  grid  boxes 
too,  which  share  the  same  ambiguous  face.  The  resulting  iso-surface  will  be  an  unique 
surface. 
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Chapter  6 


ALGORITHM  FOR  CREATION 
OF  ISO-SURFACE  PATCHES 


An  algorithm  for  connecting  the  edge  intersections  to  construct  the  polygonal 
iso-surface  patches  has  been  developed. 

6.1  The  Logic  of  the  Algorithm 

The  algorithm  treats  every  grid  box  identically.  It  starts  from  investigating  each 
of  6  faces  of  a  box.  For  each  face,  the  algorithm  finds  edge  intersections,  if  any, 
by  comparing  property  values  at  the  vertices  against  the  iso-value  Pc .  Then,  it 
establishes  the  connecting  information  among  the  points  found  on  each  face.  If  two 
points  are  found,  the  connection  is  trivial.  If  four  points  axe  found,  it  calculates  the 
centroidal  value  to  find  the  correct  way  of  connecting  points.  Those  intersection  points 
are  saved  to  an  array  with  their  coordinates,  face  identification,  edge  identification, 
and  connecting  information.  The  local  identification  system  shown  on  Figure  3.2  is 
used.  Table  6.1  shows  an  example  of  the  data  structure  of  the  array. 

After  investigating  all  six  faces  of  the  box,  the  algorithm  sorts  the  saved  edge 
intersections  to  construct  polygons.  We  have  the  number  of  edge  intersections  with 
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Table  6.1:  Data  structure  for  tbe  algorithm 

their  coordinates,  face  id’s,  edge  id’s,  and  connecting  information  on  face  level.  But, 
each  edge  intersection  is  saved  twice  from  the  previous  process  because  one  edge  is 
shared  by  two  faces.  Starting  from  any  point,  it  moves  to  the  connecting  point  on 
the  same  face.  Then  it  connects  the  duplicated  point  on  the  same  edge.  It  travels 
along  the  sides  of  an  iso-surface  patch  until  it  returns  to  the  strating  point.  If  there 
are  points  remained,  not  connected,  it  starts  from  one  of  the  remainders  and  does 
the  same  procedure  to  find  another  patch  within  the  grid  box.  Finally  the  iso-surface 
patches  are  triangulated  for  the  rendering  process.  Figure  6.1  shows  the  sorting  pro¬ 
cess,  which  travels  along  the  boundary  of  an  iso-surface  patch.  The  previous  Table  6.1 
shows  the  data  structure  after  the  sorting  process.  The  algorithm  in  pseudo-code  is 
as  follows. 
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Figure  6.1:  Traveling  along  a  polygon 

For  all  faces  of  one  box 

Find  all  edge  intersections,  if  any 
Find  connecting  information  on  the  face 
If  2  points  :  trivial 

If  4  points  :  use  function  value  of  centroid  to  get  connecting  info. 
Save  their  coordinates  with  face  id,  edge  id, 
and  face  level  connecting  information 
Sort  edge  intersections  to  construct  polygon(s) 

Triangulate  each  polygon  found 
End 

Procedure  :  Sort  edge  intersections  to  construct  polygon(s) 
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Start  from  any  point 
Find  connecting  point  on  the  same  face 
Find  duplicated  point  on  the  same  edge 
Keep  going  until  it  returns  to  the  starting  point 
If  there  are  points  remained,  not  connected,  start  one  of  them  and 
do  the  same  work  as  above 

End 

Observations  used  throughout  the  algorithm  are  summarized  as  follows  . 

•  A  box  has  6  faces,  8  vertices,  and  12  edges. 

•  One  edge  can  have  at  most  one  intersection  with  an  iso-surface.  It  is  because 
of  the  linear  interpolation  used  along  each  edge  of  the  grid  box. 

•  One  face  can  have  0,  2,  or  4  edge  intersections. 

•  One  box  can  have  at  most  12  intersections,  one  per  12  edges  of  the  box. 

•  Within  one  grid  box,  there  may  be  more  than  one  patch  of  the  iso-surface. 

•  Iso-surface  divides  the  higher  value  region  from  the  lower  value  region.  The 
surface  has  inside  and  outside  information. 

•  Polygonal  iso-surface  patches  are  non-planar  in  roost  cases. 

6.2  Calculation  of  Surface  Normals 

Normal  vectors  of  surface  patches  are  required  for  the  rendering  process  such  as 
Gouraud  shading  or  Phong  shading.  They  can  be  calculated  by  taking  cross  product 
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of  edge  vectors  of  polygonal  patches  which  share  a  vertex.  This  method  of  calculating 
normal  vectors  gives  rough  approximation.  Average  of  normal  vectors  calculated  from 
all  polygons  sharing  a  vertex  is  usually  used  for  the  final  rendering  process. 

Alternatively,  numerical  derivatives  can  be  used,  based  upon  the  fact  that  the 
directional  derivative  at  a  point  is  the  normal  vector  of  iso-surface  at  the  point.  For 
F(x,y,  z)  =  0  is  given,  VF  is  the  normal  vector  at  the  point.  The  normal  vector 
(directional  derivative)  points  toward  the  higher  value  region,  so  inside  and  outside 
information  of  the  iso-surface  is  available. 

Gradients  are  estimated  by  centred  difference  method.  The  cases  where  intervals 
are  not  same  has  been  considered.  For  z  coordinate  direction,  the  forward  difference 
is 

MXi)  - a Tj - 

The  backward  difference  is 

x  /(*<)“/(*«- 1) 

MXi)  ~  Ax* 

The  central  difference  is 

Ax,  .  /;(*,)  +  Axt  •  /;(z,) 

Azi  +  A  x/ 

Figure  6.2  shows  the  process  of  taking  the  numerical  gradient.  The  other  y,  z  com¬ 
ponents  of  VF  can  be  obtained  by  the  same  procedure. 

At  boundaries  of  the  computational  domain,  the  central  differences  cannot  be 
calculated.  Forward  differences  or  backward  differences  are  used  as  the  normal  vectors 
instead.  These  normal  vectors  are  then  normalized  against  its  own  length  to  get  the 
unit  normal  vectors. 
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Figure  6.2:  Numerical  differentiations 

6.3  Implementation  of  the  Algorithm 

Hierarchical  data  structure  is  used  for  handling  of  the  adjacency  information  of 
one  box.  See  Appendix  B.  Topological  queries  which  ask  neighboring  faces,  edges, 
and  vertices  can  be  answered  through  this  data  structure. 

Implementation  of  the  algorithm  is  done  on  a  Stellar  GS1000  series  graphics  super 
computer.  Phigs+  implementation  on  the  Stellar  Computer,  which  is  based  on  X- 
window  system,  is  used  throughout  as  the  main  graphics  package.  Fortran77  is  the 
main  programming  language.  Stellar’s  implementation  of  Phigs+  is  not  complete, 
and  its  baseline  graphics  are  X-window  system  written  in  language  C,  and  a  few  C 
procedures  are  used  too.  In  ter- language  calling  conventions  for  C  and  Fortran,  which 
is  set  by  Stellar  Computer,  are  used. 

After  the  polyhedral  approximation  of  the  iso-surface  is  constructed,  rendering  of 
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the  surface  is  handled  by  the  general  purpose  visualization  system  -  AVS  (Application 
Visualization  System)  -  supplied  by  Stellar  Computer. 
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Chapter  7 


APPLICATION  AND 
CONCLUSION 


The  developed  algorithm  can  be  used  for  &  range  of  areas  as  long  as  the  required 
data  set  can  be  supplied.  It  can  be  used  to  visualize  fluid  flows  from  the  result 
of  Computational  Fluid  Mechanics,  Meteorological  data  such  as  clouds  and  storm, 
Seismic  data  for  mineral  resource  development,  Medical  imaging  for  human  organs, 
...  etc. 

As  the  iso-surface  has  inside  -  outside  information  inheritantly,  the  developed  al¬ 
gorithm  can  be  extended  to  reconstruct  solid  object.  That  is,  topological  neighboring 
information  required  by  a  solid  modeling  system  can  be  extracted  through  application 
of  the  method  to  a  regularly  sampled  data  set.  One  good  example  is  reconstruction 
of  a  human  organ  from  a  CT  (Computed  Tomograph)  scan  data  set. 

Some  assumptions  are  made  in  the  development  of  the  algorithm  based  on  ob¬ 
servations.  They  should  be  tested  or  proved  rigorously.  Those  assumptions  are  as 
follows. 

Iso-surface  divide  hot  region  (higher  value  region)  from  cool  region  (lower  value 
region)  completely.  There  should  not  be  any  holes  on  the  iso-surface  between  above 
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two  regions,  if  continuity  of  function  in  the  domain  is  assumed.  At  one  instant  of 
time,  the  iso-surface  divide  the  two  regions  like  as  a  solid  surface. 

Can  the  algorithm  completly  covers  the  Iso-surface  ?  As  the  method  use  only 
local  information  available  on  one  grid  box,  and  does  not  consider  any  neighboring 
information  or  globality,  its  ability  to  reconstruct  the  whole  iso-surface  completely 
can  be  suspected.  Only  observation  which  enforce  the  global  connectivity  is  that  two 
neighboring  boxes  use  the  same  information  (intersection  points  and  connection  of 
them).  The  same  information  is  used  by  the  neighboring  box  for  the  shared  face.  See 
Figure  3.4 

The  developed  algorithm  is  relatively  simple  in  that  only  local  information  is  used 
for  the  construction  of  iso-surface.  This  characteristic  can  be  used  for  speeding  up 
via  parallel  processing.  The  algorithm  is  also  found  to  be  robust. 
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Appendix  A 

256  cases  of  Box-Surface 
Intersections 


All  the  256  cases  of  box-surface  intersections  are  enumerated.  They  are  grouped 
by  the  number  of  edge  intersections. 
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Appendix  A 


Tha  following  pages  show  all  the  256  cases  of  box-surface  intersection.  The  cases  are 
sorted  by  the  number  of  edge  Intersections  (IN)  resulting  in  9  different  groups. 


Explanation: 

+  indicates  P  >  Pc  at  vertex 

indicates  P  <  Pc  at  vertex 
*  indicates  edge  intersection  point 

ID  =  n=bbbbbbbb 

n  -  case  ID  number  (between  0  and  255) 
bbbbbbbb  «  corresponding  bit  pattern 
IN  number  of  edge  intersections 


Directory: 


Group 

IN 

No.  of  cases 

pages 

1 

0 

2 

A 2. 

2 

3 

16 

A.3 

3 

4 

30 

A.4  -  A.5 

4 

5 

48 

A.6  *  A.7 

5 

6 

64 

A.8  -  A.10 

6 

7 

48 

A.11  -  A.12 

7 

8 

30 

A.1 3  *  A.1 4 

8 

9 

16 

A.1 5 

9 

12 

2 

A.16 

A.1 


ID*  15=000011 1 1  IN=  4  ID*  17=00010001  IN*  4  ID*  24=00100010  IN*  4  ID*  48=00110000  IN*  4 


10*102*01100110  IN*  4  10*111*01101111  IN*  4  10-119*01110111  IN*  4  10*126*10001000  IN*  4 


10-192-11000000  IN-  4  10-204-11001100  IN-  4  10-207-11001111  IN-  4  10-221-11011101  IN-  4 


A.4 


ID-241-11110001  IN-  5  10-242-1 11 10010  IN-  5  ID-244-11110100  IN-  5  ID-241-11111090  IN- S 


A.7 


10=198=11000110  IN=  6  10=201=11001001  IN*  6  10=209=11010001  IN*  6  10=212=11010100  IN=  6 


10=215=11010111  IN*  6  10=216=11011000  IN-  6  10=219=11011011  IN-  6  10=222=11011110  IN-  6 


10-226-11100010  IN-  6  10-228-11100100  IN-  6  10-231-11100111  IN-  6  ID-232-11101000  IN-  6 


10-235-11101011  IN-  6  10-237-11101101  IN-  6  10-245-11110101  IN-  6  10-250-11111010  IN-  6 


10-109-01101101  IN-  7  10-117-01110101  IN-  7  10-121-01111001  IN-  7  10-124-01111100  IN-  7 


A.11 


A.12 


ID-105-01101001  IN-  8  ID-106-01101010  IN-  8  ID-120-01111000  IN-  8  ID-135-10000111  IN-  8 


ID-149-10010101  IN-  8  ID-150-10010110  IN-  8  ID-154-10011010  IN-  8  ID-163-10100011  IN-  8 


ID-166-10100110  IN-  8  ID-169-10101001  IN-  I  ID-170-10101010  IN-  I  10-172-10101100  IN-  • 


Appendix  B 

Hierarchical  Data  Structure 

Several  different  data  structures  are  being  used  to  represent  topological  neighbor¬ 
ing  relationship  in  the  boundary  representation  of  Solid  modeling  systems.  To  supply 
the  adjacency  information  within  a  box,  the  simple  Hierarchical  data  structure  is 
used  in  the  algorithm.  For  a  given  edge,  as  an  example,  the  algorithm  requires  the 
information  of  two  faces  which  share  the  edge,  and  two  vertices  which  are  limiting 
the  edge. 

The  Hierarchical  data  structure  treats  all  the  faces  of  the  polyhedra  as  roots.  Fur¬ 
thermore,  each  face  has  information  of  surrounding  edges.  Each  edge  has  information 
of  two  bounding  vertices  as  its  children.  The  structure  is  composed  of  2  tables,  Face- 
Edge  Table  (FET)  and  Edge- Vertex  Table  (EVT).  Figure  B.l  shows  those  relations, 
where  following  definitions  are  used. 

•  F  :  face 

•  E : edge 

•  V  :  vertex 

•  nf :  total  no.  of  faces 
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V 


F-E  table  (FED 


E-V  table  (EVP 


VI 


El 


Ere 


Figure  B.l:  Relationships  in  the  Hierarchical  Data  Structure 

•  ne  :  total  no.  of  edges 

•  nv  :  total  no.  of  vertices 

•  nee  :  no.  of  adjacent  edges,  given  E 

•  nw  :  no.  of  adjacent  vertices,  given  V 

•  nve  :  no.  of  adjacent  edges,  given  V 

•  nvf  :  do.  of  adjacent  faces,  given  V 

•  nfe  :  no.  of  adjacent  edges,  given  F 

•  nfv  :  no.  of  adjacent  vertices,  given  F 

•  nff  :  no.  of  adjacent  faces,  given  F 
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V2 


Generally,  the  data  structure  will  be  used  to  answer  the  following  9  topological  queries 


•  T1  :  Given  V,  find  all  nvv  vertices  around  V 

•  T2  :  Given  V,  find  all  nve  edges  around  V 

•  T3  :  Given  V,  find  all  nvf  faces  around  V 

•  T4  :  Given  E,  find  2  vertices  connected  by  E 

•  T5  :  Given  E,  find  all  nee  edges  around  E 

•  T6  :  Given  E,  find  2  faces  that  meet  at  E 

•  T7  :  Given  F,  find  all  nfv  vertices  around  F 

•  T8  :  Given  F,  find  all  nfe  edges  around  F 

•  T9  :  Given  F,  find  all  nff  faces  around  F 

For  a  box,  the  Hierarchical  data  structure  has  a  fixed  form.  If  we  use  the  local 

identification  system  of  Figure  3.2,  the  two  tables  have  the  form  as  Table  B.l  and 
Table  B.2. 

Another  table  is  created  to  get  the  grid  position  from  the  base  vertex  which  is  the 
grid  point  with  the  lowest  (*,  j,  k)  tuple  within  the  box.  The  table  has  the  differences 
of  grid  positions  between  each  vertex  and  the  base  vertex.  Table  B.3  shows  the 
relation. 
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■ 

Ex 

Fa 

Fa 

f4 

Fx 

1 

2 

3 

4 

Fa 

1 

9 

5 

10 

f3 

2 

11 

6 

10 

f4 

3 

12 

7 

11 

B 

4 

9 

8 

12 

H 

5 

6 

7 

8 

Table  B.l:  Face- Edge  Table 
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vx  v3 


E 

wm 

E, 

.  1 

Ex 

■ 

e4 

n 

Es 

5  6 

E 

6  7 

Er 

7  8 

E 

8  5 

Ex 

1  5 

Ex  o 

2  6 

Ex 

3  7 

E\3 

4  8 

Table  B.2:  Edge- Vertex  Table 
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■ 

Ai 

Aj 

Ak 

Vi 

0 

0 

0 

V2 

1 

D 

0 

V3 

1 

i 

V4 

0 

i 

D 

Vs 

0 

0 

i 

V* 

1 

0 

i 

V7 

1 

1 

i 

V. 

0 

1 

i 

Table  B.3:  Vertex  Difference  Table 
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